20210112 - raw track analysis notebook

A first stab at analyzing the raw data. Most important functions are placed in utils/trx_utils.py. Note that here we focus on two sets of data here: 'exp1_cam1' which is the first camera from our long timescale data and '24h_bright' which is a dataset I gathered last weekend with significantly increased brightness and a verified light-dark cycle.

import logging
from seaborn.distributions import distplot
import seaborn as sns
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from tqdm import tqdm
import pandas as pd
import h5py
import numpy as np
from pathlib import Path
import os
from datetime import datetime
import importlib

wd = "/Genomics/ayroleslab2/scott/long-timescale-behavior/analysis/"
os.chdir(wd)

import utils.trx_utils as trx_utils

data_dir = "/Genomics/ayroleslab2/scott/long-timescale-behavior/data/"
track_dir = "/Genomics/ayroleslab2/scott/long-timescale-behavior/data/tracks/"
plots_dir = "/Genomics/ayroleslab2/scott/long-timescale-behavior/analysis/plots/"

logging.basicConfig(
    format="%(asctime)s %(levelname)s: %(message)s",
    level=logging.INFO,
    datefmt="%H:%M:%S",
)
logger = logging.getLogger("analysis_logger")

Metadata

Build up a dict containing all of our experimental metadata so we can easily access it later.

exp1_cam1_h5s = [
    "exp2_cam1_0through23.tracked.analysis.h5",
    "exp2_cam1_24through47.tracked.analysis.h5",
    "exp2_cam1_48through71.tracked.analysis.h5",
    "exp2_cam1_72through95.tracked.analysis.h5",
    "exp2_cam1_96through119.tracked.analysis.h5",
]
exp1_cam1_h5s = [track_dir + filename for filename in exp1_cam1_h5s]

bright_h5s = ["24h_bright_0through23_0through29.tracked.analysis.h5"]
bright_dir = "/Genomics/ayroleslab2/scott/long-timescale-behavior/tmp/24h_bright/"
bright_h5s = [bright_dir + filename for filename in bright_h5s]


FMT = "%w-%H:%M:%S"

# Build with dict for compatibility with JSON
expmt_dict = {
    "exp1_cam1": {
        "h5s": exp1_cam1_h5s,
        "video": "exp1_cam1.mkv",
        "frame_rate": 100,
        "start_time": datetime.strptime("0-22:33:00", FMT),
        "camera": "1",
        "experiment": "1",
        "video_path": "/Genomics/ayroleslab2/scott/long-timescale-behavior/data/exp1/exp5_202109014_2233/Camera1/exp.mkv",
        "px_mm": 28.25,
    },
    "24h_bright": {
        "h5s": bright_h5s,
        "video": "/Genomics/ayroleslab2/scott/long-timescale-behavior/tmp/24h_bright/24h_bright.mkv",
        "frame_rate": 99.96,
        "start_time": datetime.strptime("0-12:00:00", FMT),
        "camera": "1",
        "experiment": "24h_bright",
        "video_path": "/Genomics/ayroleslab2/scott/long-timescale-behavior/tmp/24h_bright/24h_bright.mkv",
        "px_mm": 28.25,
    },
}
px_mm = expmt_dict["exp1_cam1"]["px_mm"]

# Create a dict for placing tracks. We use dicts to separate data and metadata without needing to move into classes.
tracks_dict_raw = {}

Load h5 traces and match by quadrant

For reference, fly node locations are in the form (time, node, coord, fly_idx). The frequency and assignments are logged here but are not necessary for understanding the subsequent analysis. We expect 120*60*60*100 = 43,200,000 frames for first 120h of exp1_cam1 and 24*60*60*99.96 = 8,636,544 frames for 24h_bright. The missing values here are just locations where the thorax was not identified.

# For each experiment, load the tracks, sort, merge, and add to track_dict_raw
for key in expmt_dict:
    expmt_name = str(key)
    logger.info(f"Loading {expmt_name}")
    expmt = expmt_dict[key]

    with h5py.File(expmt["h5s"][0], "r") as f:
        logger.info(expmt["h5s"][0])
        dset_names = list(f.keys())
        # Note the assignment of node_names here!
        node_names = [n.decode() for n in f["node_names"][:]]
        expmt_dict[key]["node_names"] = node_names
        locations = f["tracks"][:].T

        locations[:, :, 1, :] = -locations[:, :, 1, :]
        assignment_indices, locations, freq = trx_utils.hist_sort(
            locations, ctr_idx=node_names.index("thorax"), ymin=-1536, ymax=0
        )
        locations[:, :, 1, :] = -locations[:, :, 1, :]

    if len(expmt["h5s"]) > 1:
        for filename in tqdm(expmt["h5s"][1:]):
            with h5py.File(filename, "r") as f:
                temp_locations = f["tracks"][:].T
                temp_locations[:, :, 1, :] = -temp_locations[:, :, 1, :]
                temp_assignment_indices, temp_locations, freq = trx_utils.hist_sort(
                    temp_locations,
                    ctr_idx=node_names.index("thorax"),
                    ymin=-1536,
                    ymax=0,
                )
                temp_locations[:, :, 1, :] = -temp_locations[:, :, 1, :]

                logger.info(filename)
                logger.info(freq)

                locations = np.concatenate((locations, temp_locations), axis=0)

    # Final assignment
    locations[:, :, 1, :] = -locations[:, :, 1, :]
    assignment_indices, locations, freq = trx_utils.hist_sort(
        locations, ctr_idx=node_names.index("thorax"), ymin=-1536, ymax=0
    )
    locations[:, :, 1, :] = -locations[:, :, 1, :]
    logger.info(f"Experiment: {str(expmt)}")
    logger.info(f"Final frequencies: {freq}")
    logger.info(f"Final assignments: {assignment_indices}")

    tracks_dict_raw[expmt_name] = locations  # [0:1000,:,:,:]
    expmt_dict[expmt_name]["assignments"] = assignment_indices
    expmt_dict[expmt_name]["freq"] = freq

23:00:40 INFO: Loading exp1_cam1
23:00:40 INFO: /Genomics/ayroleslab2/scott/long-timescale-behavior/data/tracks/exp2_cam1_0through23.tracked.analysis.h5
  0%|          | 0/4 [00:00<?, ?it/s]23:01:35 INFO: /Genomics/ayroleslab2/scott/long-timescale-behavior/data/tracks/exp2_cam1_24through47.tracked.analysis.h5
23:01:35 INFO: [array([[      0., 8639891.],
       [      0.,       0.]]), array([[      0.,       0.],
       [      0., 8638655.]]), array([[8639782.,       0.],
       [      0.,       0.]]), array([[      0.,       0.],
       [8639966.,       0.]])]
 25%|██▌       | 1/4 [00:29<01:29, 29.86s/it]23:02:02 INFO: /Genomics/ayroleslab2/scott/long-timescale-behavior/data/tracks/exp2_cam1_48through71.tracked.analysis.h5
23:02:02 INFO: [array([[      0.,       0.],
       [      0., 8638097.]]), array([[      0., 8639751.],
       [      0.,       0.]]), array([[      0.,       0.],
       [8639876.,       0.]]), array([[8639751.,       0.],
       [      0.,       0.]])]
 50%|█████     | 2/4 [00:59<00:59, 29.56s/it]23:02:33 INFO: /Genomics/ayroleslab2/scott/long-timescale-behavior/data/tracks/exp2_cam1_72through95.tracked.analysis.h5
23:02:33 INFO: [array([[      0.,       0.],
       [      0., 8639624.]]), array([[      0., 8639474.],
       [      0.,       0.]]), array([[      0.,       0.],
       [8639951.,       0.]]), array([[8639662.,       0.],
       [      0.,       0.]])]
 75%|███████▌  | 3/4 [01:32<00:31, 31.41s/it]23:03:05 INFO: /Genomics/ayroleslab2/scott/long-timescale-behavior/data/tracks/exp2_cam1_96through119.tracked.analysis.h5
23:03:05 INFO: [array([[      0.,       0.],
       [      0., 8639434.]]), array([[      0., 8639595.],
       [      0.,       0.]]), array([[      0.,       0.],
       [8639888.,       0.]]), array([[8639866.,       0.],
       [      0.,       0.]])]
100%|██████████| 4/4 [02:05<00:00, 31.37s/it]
23:03:35 INFO: Experiment: {'h5s': ['/Genomics/ayroleslab2/scott/long-timescale-behavior/data/tracks/exp2_cam1_0through23.tracked.analysis.h5', '/Genomics/ayroleslab2/scott/long-timescale-behavior/data/tracks/exp2_cam1_24through47.tracked.analysis.h5', '/Genomics/ayroleslab2/scott/long-timescale-behavior/data/tracks/exp2_cam1_48through71.tracked.analysis.h5', '/Genomics/ayroleslab2/scott/long-timescale-behavior/data/tracks/exp2_cam1_72through95.tracked.analysis.h5', '/Genomics/ayroleslab2/scott/long-timescale-behavior/data/tracks/exp2_cam1_96through119.tracked.analysis.h5'], 'video': 'exp1_cam1.mkv', 'frame_rate': 100, 'start_time': datetime.datetime(1900, 1, 1, 22, 33), 'camera': '1', 'experiment': '1', 'video_path': '/Genomics/ayroleslab2/scott/long-timescale-behavior/data/exp1/exp5_202109014_2233/Camera1/exp.mkv', 'px_mm': 28.25, 'node_names': ['head', 'eyeL', 'eyeR', 'thorax', 'abdomen', 'forelegL', 'forelegR', 'midlegL', 'midlegR', 'hindlegL', 'hindlegR', 'wingL', 'wingR', 'proboscis']}
23:03:35 INFO: Final frequencies: [array([[43199015.,        0.],
       [       0.,        0.]]), array([[       0., 43198678.],
       [       0.,        0.]]), array([[       0.,        0.],
       [43199650.,        0.]]), array([[       0.,        0.],
       [       0., 43195769.]])]
23:03:35 INFO: Final assignments: [0 1 2 3]
23:03:35 INFO: Loading 24h_bright
23:03:35 INFO: /Genomics/ayroleslab2/scott/long-timescale-behavior/tmp/24h_bright/24h_bright_0through23_0through29.tracked.analysis.h5
23:04:12 INFO: Experiment: {'h5s': ['/Genomics/ayroleslab2/scott/long-timescale-behavior/tmp/24h_bright/24h_bright_0through23_0through29.tracked.analysis.h5'], 'video': '/Genomics/ayroleslab2/scott/long-timescale-behavior/tmp/24h_bright/24h_bright.mkv', 'frame_rate': 99.96, 'start_time': datetime.datetime(1900, 1, 1, 12, 0), 'camera': '1', 'experiment': '24h_bright', 'video_path': '/Genomics/ayroleslab2/scott/long-timescale-behavior/tmp/24h_bright/24h_bright.mkv', 'px_mm': 28.25, 'node_names': ['head', 'eyeL', 'eyeR', 'thorax', 'abdomen', 'forelegL', 'forelegR', 'midlegL', 'midlegR', 'hindlegL', 'hindlegR', 'wingL', 'wingR', 'proboscis']}
23:04:12 INFO: Final frequencies: [array([[8636632.,       0.],
       [      0.,       0.]]), array([[      0., 8636601.],
       [      0.,       0.]]), array([[      0.,       0.],
       [8625652.,       0.]]), array([[      0.,       0.],
       [      0., 8636473.]])]
23:04:12 INFO: Final assignments: [0 1 2 3]

Process raw tracks

Processing the raw tracks as needed and generating corresponding velocities. Note that the velocities use np.gradient which uses the central difference for interior points and the forward difference for the first and last points.

importlib.reload(trx_utils)
tracks_dict = {}
velocities_dict = {}

for key in expmt_dict:
    expmt = expmt_dict[key]
    fly_node_locations = tracks_dict_raw[key].copy()
    fly_node_locations = trx_utils.fill_missing_np(fly_node_locations)

    # fly_node_locations = trx_utils.smooth_median(fly_node_locations, window=5)
    # fly_node_locations = trx_utils.smooth_gaussian(fly_node_locations, window=5)

    fly_node_velocities = (
        trx_utils.instance_node_velocities(
            fly_node_locations, 0, fly_node_locations.shape[0]
        )
        * (1 / px_mm)
        * expmt["frame_rate"]
    )
    tracks_dict[key] = fly_node_locations  # _filled
    velocities_dict[key] = fly_node_velocities

100%|██████████| 112/112 [06:03<00:00,  3.24s/it]
100%|██████████| 14/14 [00:41<00:00,  2.94s/it]
100%|██████████| 14/14 [00:37<00:00,  2.68s/it]
100%|██████████| 14/14 [00:35<00:00,  2.54s/it]
100%|██████████| 14/14 [00:40<00:00,  2.86s/it]
100%|██████████| 112/112 [01:02<00:00,  1.79it/s]
100%|██████████| 14/14 [00:07<00:00,  1.97it/s]
100%|██████████| 14/14 [00:06<00:00,  2.15it/s]
100%|██████████| 14/14 [00:06<00:00,  2.14it/s]
100%|██████████| 14/14 [00:06<00:00,  2.15it/s]

Save

Save the JSON and h5s if needed. We'll avoid using high compression options to keep our time cost low.

import json
json.dump(expmt_dict, open("expmt_dict.json", "w"), default=str)

for key in tqdm(expmt_dict):
    data_file = h5py.File(data_dir + f"{key}_fly_node_locations.h5", "w")
    data_file.create_dataset(
        "tracks", data=tracks_dict[key]
    )  # , compression='lzf')#'gzip', compression_opts=9)
    data_file.close()

    data_file = h5py.File(data_dir + f"{key}_fly_node_velocities.h5", "w")
    data_file.create_dataset(
        "velocities", data=velocities_dict[key]
    )  # , compression='lzf')#'gzip', compression_opts=9)
    data_file.close()

100%|██████████| 2/2 [01:35<00:00, 47.85s/it]

Load in processed tracks

This saves a good chunk of time by skipping the first section of the notebook.

import json

expmt_dict = json.load(open("expmt_dict.json", "r"))
node_names = expmt_dict["exp1_cam1"]["node_names"]
px_mm = expmt_dict["exp1_cam1"]["px_mm"]
tracks_dict = {}
velocities_dict = {}
for key in tqdm(expmt_dict):
    with h5py.File(data_dir + f"/{key}_fly_node_locations.h5", "r") as f:
        tracks_dict[key] = f["tracks"][:]
    with h5py.File(data_dir + f"/{key}_fly_node_velocities.h5", "r") as f:
        velocities_dict[key] = f["velocities"][:]

100%|██████████| 2/2 [02:43<00:00, 81.82s/it] 

Time series check

Note that we need to be super careful about the axes as they vary between plots. We start off by plotting a minute start from the beggining of the experiment where all four individuals are not moving. It appears that at least some portion of our data is fairly quantized For example, when looking at the x position of the first track (exp1_cam1, fly0) we can see the repeated oscillation between ~172.1 and ~172.8. However, if we take a look at the x coordinate of the second track (exp1_cam1, fly1) we can see that the data in a much tighter band between ~522.9 and 523.1.

importlib.reload(trx_utils)
expmt_name = "exp1_cam1"
frame_start = 0
frame_end = int(frame_start + 60 * 100)
vels = velocities_dict[expmt_name]

# trx_utils.plot_timeseries(tracks_dict[expmt_name],vels,node_idx=node_names.index("thorax"),frame_start=frame_start,frame_end=frame_end, title = expmt_name,output_name = f'{expmt_name}_thorax_timeseries.png',path='')
trx_utils.plot_timeseries(
    tracks_dict[expmt_name],
    [0, 1, 2, 3],
    vels=None,
    node_idx=node_names.index("thorax"),
    frame_start=frame_start,
    frame_end=frame_end,
    title=expmt_name,
    output_name=f"{expmt_name}_thorax_timeseries.png",
    path="",
)

Now lets take a look at when some individuals are moving. Here we take a look at another 60 second window when three individuals are moving.

import time
expmt_name = "exp1_cam1"
frame_start = 20000 * 100
frame_end = int(frame_start + 60 * 100)
trx_utils.plot_timeseries(
    tracks_dict[expmt_name],
    [0, 1, 2, 3],
    vels=None,
    node_idx=node_names.index("thorax"),
    frame_start=frame_start,
    frame_end=frame_end,
    title=expmt_name,
    output_name=f'{time.strftime("%Y%m%d_%H%M%S")}_{frame_start}to{frame_end}_{expmt_name}_thorax_timeseries.png',
    path="",
)

Now lets take a closer look at the x coordinate of the first track (exp1_cam1, fly0) and see what it looks like when the fly is moving. Here we zoom in to a 1/10 sec region. Notice the much smoother transitions.

importlib.reload(trx_utils)
expmt_name = "exp1_cam1"
frame_start = 20055 * 100
frame_end = int(frame_start + 0.05 * 100)
trx_utils.plot_timeseries(
    tracks_dict[expmt_name],
    [0],
    vels=None,
    node_idx=node_names.index("thorax"),
    frame_start=frame_start,
    frame_end=frame_end,
    title=expmt_name,
    output_name=f'{time.strftime("%Y%m%d_%H%M%S")}_{frame_start}to{frame_end}_{expmt_name}_thorax_timeseries_zoomed.png',
    path="",
)

What happens if we dig into the highest velocity locations?

The tracks below show the timeseries at the max velocity.

expmt_name = "exp1_cam1"
fly_idx = 0
width = 100
node = node_names.index("thorax")
vel_max = np.argmax(velocities_dict[expmt_name][:, node, fly_idx])
importlib.reload(trx_utils)
logger.info(f'"Plotting fly {fly_idx}, node {node_names[node]}, max vel {velocities_dict[expmt_name][vel_max, node, fly_idx]} at frame {vel_max}')
frame_start = vel_max - width
frame_end = int(frame_start +  2*width -1)
trx_utils.plot_timeseries(
    tracks_dict[expmt_name],
    [fly_idx],
    vels=None,
    node_idx=node,
    frame_start=frame_start,
    frame_end=frame_end,
    title=expmt_name,
    output_name=f'{time.strftime("%Y%m%d_%H%M%S")}_{frame_start}to{frame_end}_{expmt_name}_thorax_timeseries_zoomed.png',
    path="",
)

width=20
logger.info(f'"Plotting fly {fly_idx}, node {node_names[node]}, max vel {velocities_dict[expmt_name][vel_max, node, fly_idx]} at frame {vel_max}')
frame_start = vel_max - width
frame_end = int(frame_start +  2*width-1)
trx_utils.plot_timeseries(
    tracks_dict[expmt_name],
    [fly_idx],
    vels=None,
    node_idx=node,
    frame_start=frame_start,
    frame_end=frame_end,
    title=expmt_name,
    output_name=f'{time.strftime("%Y%m%d_%H%M%S")}_{frame_start}to{frame_end}_{expmt_name}_thorax_timeseries_zoomed.png',
    path="",
)

23:20:04 INFO: "Plotting fly 0, node thorax, max vel 414.882712974806 at frame 9576770

23:20:04 INFO: "Plotting fly 0, node thorax, max vel 414.882712974806 at frame 9576770

Speed distributions

Here we plot the speed distributions -- central difference except on edges calculated by np.gradient -- for exp1_cam1 and 24h_bright. One very concerning aspect is that we see drastically more movement in the 24h_bright data than in the exp1_cam1 data. We can also see the tendency to move in pixel steps just as we would expect from the quantization described above..

import palettable

for expmt_id in np.arange(len(list(expmt_dict.keys()))):
    expmt_name = list(expmt_dict.keys())[expmt_id]
    for fly_id in np.arange(tracks_dict[expmt_name].shape[3]):
        fig, ax = plt.subplots()
        ax.hist(
            velocities_dict[expmt_name][:, node_names.index("thorax"), fly_id],
            bins=np.arange(0, 20, 0.25),
            color=palettable.wesanderson.Zissou_5.mpl_colors[expmt_id],
        )
        ax.set_title(f"Histogram of thorax velocities for {expmt_name} - fly {fly_id}")
        ax.set_xlabel("Velocity (mm/s)")
        fig.savefig(
            f'{time.strftime("%Y%m%d_%H%M")}_{expmt_name}_fly{fly_id}_thorax_velocities.png'
        )

Video plotting

Good for documentation and debugging but ignored until we're happy with timeseries.

%%capture
# Plot traces on individuals/video as needed

import time
importlib.reload(trx_utils)
expmt_name = 'exp1_cam1'
frame_start = 20025*100
frame_end = int(frame_start + 60*100)

trx_utils.plot_trx(tracks_dict[expmt_name],expmt_dict[expmt_name]["video_path"],frame_start,frame_end,output_path=f'{time.strftime("%Y%m%d_%H%M%S")}_{frame_start}to{frame_end}_{expmt_name}_trx')

from IPython.display import Video
import os

cwd = os.getcwd()
Video(f"{cwd}/plots/20220112_232029_2002500to2008500_exp1_cam1_trx_fly_node_locations.mp4",embed=True)

# A simple example of the plot_ego function

importlib.reload(trx_utils)
expmt_name = "24h_bright"
frame_start = int(27145 * 100)
frame_end = int(27145 * 100 + 600 * 100)

for fly_id in range(tracks_dict[expmt_name].shape[3]):
    egocentric_node_locations, egocentric_angles = trx_utils.normalize_to_egocentric(
        tracks_dict[expmt_name][:, :, :, fly_id],
        ctr_ind=node_names.index("thorax"),
        fwd_ind=node_names.index("head"),
        return_angles=True,
    )
    # egocentric_angles=None
    trx_utils.plot_ego(
        tracks_dict[expmt_name],
        expmt_dict[expmt_name]["video_path"],
        egocentric_angles,
        [fly_id],
        node_names.index("thorax"),
        frame_start,
        frame_end,
        output_path=f'{time.strftime("%Y%m%d_%H%M%S")}_{frame_start}to{frame_end}_{expmt_name}_fly{fly_id}_raw_ego.mp4',
    )
    egocentric_velocities = (
        trx_utils.instance_node_velocities(
            egocentric_node_locations, 0, egocentric_node_locations.shape[0]
        )
        * (1 / px_mm)
        * expmt["frame_rate"]
    )

Circadian rhythm

Good for documentation -- to be ignored for now.

import scipy.stats
from datetime import datetime

ToD = {}
for expmt in expmt_dict:
    expmt_data = expmt_dict[expmt]
    FMT = "%w-%H:%M:%S"
    start_day = datetime.strptime(
        "1900-01-01 00:00:00", "%Y-%m-%d %H:%M:%S"
    )  # for example
    try:
        time = datetime.strptime(expmt_data["start_time"], "%Y-%m-%d %H:%M:%S")
    except:
        time = expmt_data["start_time"]
    expmt_data["start_time"] = time
    difference = start_day - time
    frame_rate = expmt_data["frame_rate"]
    shift = int(difference.seconds * frame_rate)
    frame_idx = np.arange(tracks_dict[expmt].shape[0]) - shift
    expmt_tod = (frame_idx % int(24 * 60 * 60 * frame_rate)) / (
        1 * 60 * 60 * frame_rate
    )
    ToD[expmt] = expmt_tod

day_start = 8
day_end = 20
day_dict = {}
for expmt in expmt_dict:
    day_dict[expmt] = (ToD[expmt] > day_start) & (ToD[expmt] < day_end)

plt.rcParams["patch.linewidth"] = 0
plt.rcParams["patch.edgecolor"] = "none"
plt.rcParams["figure.figsize"] = (9, 3)
plt.rcParams["figure.dpi"] = 300

for expmt in expmt_dict:
    vels = velocities_dict[expmt].copy()
    day = (ToD[expmt] > day_start) & (ToD[expmt] < day_end)
    for fly_idx in range(vels.shape[2]):
        fly_thorax_vel = vels[:, node_names.index("thorax"), fly_idx]
        # fly_thorax_vel[fly_thorax_vel < 3*(1/28.25)*100] = 0
        binned = scipy.stats.binned_statistic(
            day.astype(int),
            fly_thorax_vel,
            statistic="mean",
            bins=[0, 0.5, 1],
            range=None,
        )
        sns.barplot(x=["Night", "Day"], y=binned.statistic)
        plt.title(f"{expmt} - Fly {fly_idx} - Mean thorax velocity by day/night")
        plt.show()

        sns.boxplot(x=day.astype(int), y=fly_thorax_vel, showfliers=False)
        plt.title(f"{expmt} - Fly {fly_idx} - Mean thorax velocity by day/night")
        plt.show()
        logger.info(f"{expmt} {fly_idx}")
        segments = 2 * 24

        binned = scipy.stats.binned_statistic(
            ToD[expmt], fly_thorax_vel, statistic="mean", bins=segments, range=None
        )
        custom_params = {"axes.spines.right": False, "axes.spines.top": False}
        sns.set_theme(style="ticks", rc=custom_params)
        fig, ax = plt.subplots()
        sns.barplot(
            ax=ax,
            x=np.arange(segments),
            y=binned.statistic,
            color=palettable.wesanderson.GrandBudapest4_5.mpl_colors[0],
        )  # ,alpha=1,width=1)
        plt.xticks([i * (segments // 24) for i in [0, 8, 20, 24]], [0, 8, 20, 24])
        plt.tight_layout(pad=2)
        # ax.set_yscale('log')
        def change_width(ax, new_value):
            for patch in ax.patches:
                current_width = patch.get_width()
                diff = current_width - new_value

                # we change the bar width
                patch.set_width(new_value)

                # we recenter the bar
                patch.set_x(patch.get_x() + diff * 0.5)

        plt.xlabel("Hour of the day")
        plt.ylabel("Mean thorax velocity (mm/s)")
        plt.title("Mean thorax velocity by hour of the day")
        change_width(ax, 0.95)
        plt.savefig(f"{expmt}_fly{fly_idx}_thorax_velocity_by_hour_of_day.png")
        plt.show()

Periodogram

Direct scipy periodogram example. I tested this on a

import scipy

fly_id = 3
expmt_name = "exp1_cam1"

vels = velocities_dict[expmt_name][:, node_names.index("thorax"), fly_id]
vels_len = vels.shape[0]

# Note that boxcar here is the same the default
f, Pxx = scipy.signal.periodogram(
    vels,
    fs=100,
    window=scipy.signal.get_window(("boxcar"), vels_len),
    scaling="density",
)
plt.semilogy(f, Pxx)
plt.title(f"{expmt_name} fly {fly_id} - scipy.periodogram")
plt.xlabel("frequency [Hz]")
plt.ylabel("PSD [V**2/Hz]")
plt.show()


# We lets check on 24h_bright as well
fly_id = 3
expmt_name = "24h_bright"

vels = velocities_dict[expmt_name][:, node_names.index("thorax"), fly_id]
vels_len = vels.shape[0]

f, Pxx = scipy.signal.periodogram(
    vels,
    fs=100,
    window=scipy.signal.get_window(("boxcar"), vels_len),
    scaling="density",
)
plt.semilogy(f, Pxx)
plt.title(f"{expmt_name} fly {fly_id} - scipy.periodogram")
plt.xlabel("frequency [Hz]")
plt.ylabel("PSD [V**2/Hz]")
plt.show()